//========================================================================================
// Mer.hh
// ---------------------------------------------------------------------------------------
//
//----------------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/include/Mer.hh#11 $
//----------------------------------------------------------------------------------------
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#ifndef _MER
#define _MER

#include "CompBioTools.hh"
#include "Sequence3Bit.hh"

#define A (0ULL)
#define C (1ULL)
#define G (2ULL)
#define T (3ULL)
#define Akjd (0ULL)
#define Ckjd (1ULL)
#define Gkjd (3ULL)
#define Tkjd (2ULL)

// want 64-bit Words 
typedef unsigned long long Word;

// want 32-bit Indexes. 
// On Itanium2, long is 64-bits, so need to define it as an uint. 
#ifdef ICC 
typedef unsigned int Index;
#else
typedef unsigned long Index;
#endif

//----------------------------------------------------------------------------------------
// * Mer
//----------------------------------------------------------------------------------------
// Class to represent an n-mer up to 32bp with a 2-bit/bp encoding. This class is primarialy
// used for hashing n-mers where keeping the total space to a minimum is important.  Since 
// systems we use have 32 bit ints, Mer uses two ints to represent a 32bp sequence.  For 
// For efficiency at I/O, the encoding is masked directly from the ASCII:
//
//      Hex         Dec            Binary   Encoding Bits
// a = 61 (hex)    97 (dec)      0110 0001   0110 0(00)1   a -> 0
// c = 63 (hex)    99 (dec)      0110 0011   0110 0(01)1   c -> 1
// g = 67 (hex)    103 (dec)     0110 0111   0110 0(11)1   g -> 3
// t = 74 (hex)    116 (dec)     0111 0100   0111 0(10)0   t -> 2

// A = 41 (hex)    65 (dec)      0100 0001   0100 0(00)1   A -> 0
// C = 43 (hex)    67 (dec)      0100 0011   0100 0(01)1   C -> 1
// G = 47 (hex)    71 (dec)      0100 0111   0100 0(11)1   G -> 3 
// T = 54 (hex)    84 (dec)      0101 0100   0101 0(10)0   T -> 2 

// So, bits 0,2 effectively code for a,c,t,g


// 3-bit encoding to include N/X.  


#define BITSPERBASE (2UL)

// When we use 3-bit bases, would like encoding for n/x to be a 
// superset of the acgtACGT encoding.   Bit 3 is set for n/x where it 
// is not set for acgtACGT, so bit 3 is diagnostic and bits 1,2,3 can 
// be used as the 3-bit encoding. 

// n = 6e (hex)    110 (dec)     0110 1110   0110 (111)0   n -> 7
// x = 78 (hex)    120 (dec)     0111 1000   0111 (100)0   x -> 4
//
// N = 4e (hex)    78 (dec)      0100 1110   0100 (111)0   N -> 7
// X = 58 (hex)    88 (dec)      0101 1000   0101 (100)0   X -> 4
//
// 
class Mer{
 public:

  // For translation to Paul's encoding. 
  int BASESUNUSED; // Default is 32-mer. 
  int BITSUNUSED;

  static int mMerSize;
  
  unsigned int mUpper32;
  unsigned int mLower32;

  Mer(){};

  Mer(string s);

  Mer(Sequence3Bit &Seq,int merStart);
  Mer(Sequence3Bit &Seq,int merStart,bool &ContainsNX);

  Mer(Sequence &Seq);
  Mer(Sequence &Seq,int merStart);
  Mer(Sequence &Seq,int merStart,bool &bContainsNX);

  
  // Methods for {A,C,G,T} = bit1,bit2 of ASCII
  void ReverseComplement(Mer &fwdMer);
  void ReverseComplement();
  inline unsigned int ComplementBits(unsigned int bits);

  inline unsigned int Base2Bit(char b);
  void ToSequence(Sequence &S);
  inline char Bit2Base(unsigned int x);

  // Return the 64-bit equivalent. 
  /* Arbitrary encoding version. 
     unsigned long long ToWord() const {
     unsigned long long temp = mUpper32;
     temp = temp << 32;
     temp = temp | mLower32;
     return(temp);
     }
  */

  void SetMerSize(int merSize){
    BASESUNUSED = 32 - merSize; 
    BITSUNUSED=(BASESUNUSED*BITSPERBASE);
    mMerSize = merSize;
  }
 
  // Return the 64-bit equivalent in Paul's encoding. 
  unsigned long long ToWord() const {

    const unsigned char transtable = // 0xb4; // 10110100 (11 -> 10, 10 -> 11, etc)
      (A << (Akjd * 2)) | (C << (Ckjd * 2)) | (G << (Gkjd * 2)) | (T << (Tkjd * 2));
    // printf("%x\n", transtable);
    Index upper32,lower32;
    upper32 = mUpper32;
    lower32 = mLower32;
		
    // cout << "Upper32: " << hex << upper32 << endl;
    // cout << "Lower32: " << lower32 << dec << endl;
    Word w = 0;
		
    Index phhup = 0;
    Index phhlo = 0;
    Index revup = 0;
    Index revlo = 0;
    char kjdbits, phhbits;
		
    // OK, this is crazy.  My Lower32 bits will become
    // the upper 32 bits for the hash, and as I shift my 
    // bits out to the right, they are added to the left. 
		
    // Shift in the lower 32 bits (2 bits at a time). 
    // Set up size for divisions

#define BITSPERINDEX (sizeof(Index) * 8)
		
    for(int i = 0;i< 16;i++){
      kjdbits = (lower32 & 0x00000003) << 1; // grab the lower two bits, * 2
      lower32 = lower32>>2; // Shift them out. 
			
      phhbits = (transtable >> kjdbits) & 0x3;
      phhup <<= BITSPERBASE;
      phhup |= phhbits;
			
      revlo >>= BITSPERBASE;
      revlo |= ((Index) phhbits) << (BITSPERINDEX - BITSPERBASE);
			
      // Shift in the upper 32 bits (2 bits at a time). 
			
      kjdbits = (upper32 & 0x00000003) << 1; // grab the lower two bits, * 2. 
      upper32 = upper32>>2; // Shift them out. 
			
      phhbits = (transtable >> kjdbits) & 0x3;
      phhlo <<= BITSPERBASE;
      phhlo |= phhbits;
			
      revup >>= BITSPERBASE;
      revup |= ((Index) phhbits) << (BITSPERINDEX - BITSPERBASE);
    }
    w = (((Word) phhup) << 32) | phhlo;
    Word r = ~((((Word) revup) << 32) | revlo);
    if (BASESUNUSED) {
      w >>= BITSUNUSED;
      r <<= BITSUNUSED;
      r >>= BITSUNUSED;
    }
    // cout<<"w= "<< hex << w << " " << Mer2CString(w) <<endl;
    // cout<<"r= "<< r << " " << Mer2CString(r) <<endl;
    w = (w <= r)? w : r;
    // cout<<"min= "<< w << dec << " " << Mer2CString(w) <<endl;
    return(w);
  }
	

  void Print(ostream &fout){
    Sequence S;
    ToSequence(S);
    fout<<S;
  }

  bool operator==(const Mer &m) const {
    return((m.mUpper32 == mUpper32) && (m.mLower32 == mLower32));
  }


  bool InSubset(unsigned int Subset);


  /* 
     bool NextMer(ifstream &fin){
     // Erase any existing sequence, just in case.  
     clear();
  
     // Must be first line  
     if (buffer[0] == '1'){
     fin.getline(buffer,bufSize);
     }
  
     // Extract name and description from line.  
     string nextDesc = buffer;
     int offset = 1; // To skip >
     // Everything up to first white space is name, after that is description. 
     mName = GetNextWSDelimitedString(nextDesc,offset,offset);  
     // It's possible there is no description, so check:
     if ((offset >= 0) && (offset < nextDesc.size())){
     mDescription = nextDesc.substr(offset,nextDesc.size()-offset);
     }else{
     mDescription = "";
     }
  
     // Get the sequence until another > is encountered. 
     do{
     fin.getline(buffer,bufSize);
     string buf = buffer;
     int bytesRead = fin.gcount();
     int bsize = strlen(buffer);
     // This means that we have a null before the newline, so I abort this
     // probably corrupt read. 
     if ((bsize < (bytesRead-1)) && (bsize !=0) && (bytesRead != 0)) return(false);
	 
     if (buffer[0] != '>'){
     // Sadly, many FASTA files have empty lines or whitespace at the 
     // end of lines.  This is effectively filtered here by adding a null
     // terminator in place of all whitespace. 
		
     register int sEnd = bsize-1;
     while((buffer[sEnd] == '\t') || 
     (buffer[sEnd] == ' ') ||
     (buffer[sEnd]=='\n') && (sEnd > 0)){
     cout<<"Consuming white space: "<<hex<<(int)buffer[sEnd]<<" sEnd="<<sEnd<<dec<<endl;
     buffer[sEnd--] = '\0';
     }

     // Copy the characters of the line to the sequence vector. 
		
     ConcatenateCString(buffer);			            
     }
     }while((!fin.eof()) && (buffer[0] != '>'));
     return(true);
     }
  */

};


struct eqMer{
  bool operator()(const Mer a,const Mer b) const {
    // The mers match in the same sense. 
    if ((a.mUpper32 == b.mUpper32) && (a.mLower32 == b.mLower32)){
      return(true);
    }else return(false);
  }
};

// An equal function that does not distinguish between fwd and 
// revcomp versions of a Mer. 
// 
struct eqMerSymmetric{

  bool operator()(const Mer a,const Mer b) const {
    // The mers match in the same sense. 
    if ((a.mUpper32 == b.mUpper32) && (a.mLower32 == b.mLower32)){
      return(true);
    }else{
      // see if the reverse complement matches. 
      Mer rcb = b;
      rcb.ReverseComplement();
      if ((a.mUpper32 == rcb.mUpper32) && (a.mLower32 == rcb.mLower32)){
	return(true);
      }else{
	return(false);
      }
    }
  }
};


#endif
